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Abstract 

We present a rigorous analysis of the role and uses of the adiabatic invariant in the 
Mixmaster dynamical system. We propose a new invariant for the global dynamics 
which in some respects has an improved behaviour over the commonly used one. We 
illustrate its behaviour in a number of numerical results. We also present a new 
formulation of the dynamics via Catastrophe Theory. We find that the change from 
one era to the next corresponds to a fold catastrophe, during the Kasner shifts the 
potential is an Implicit Function Form whereas, as the anisotropy dissipates, the 
Mixmaster potential must become a Morse 0-saddle. We compare and contrast our 
results to many known works on the Mixmaster problem and indicate how extensions 
could be achieved. Further exploitation of this formulation may lead to a clearer 
understanding of the global Mixmaster dynamics. 



I Introduction 

During the past eight years or so interest in the Mixmaster universe has increased 
dramatically. There are at least two reasons for this noticeable increase. Firstly, 
there was already a considerable amount of background work concerning the basic 
dynamical issues of the model [@, [|, f|, |5|] and it is probably fair to say that this model 
(diagonal Bianchi IX) was studied more than any other homogeneous cosmology 
(with a probable exception of inflationary issues). The picture drawn from that boby 
of work was already quite rich (the actual system was also integrated numerically 
by Zardecki || but some of his results later criticized as incompatible with those of 
BKL in 0, §]) and allowed for further generalizations to be considered. 

Examples of this kind of generalized problems included Kaluza-Klein extensions 

and the search for complicated Mixmaster behaviour in other theories of gravitation 

TDJ. All these generalizations were motivated mainly by an idea and results first 



obtained by J.D. Barrow 0] almost another decade earlier, that the well-known BKL- 
Misner oscillatory, Mixmaster behaviour should in fact be viewed as an example of 
the manifestation of chaotic (erratic, unpredictable) structures in general relativity. 
It was therefore natural to examine how common such a behaviour could be. (It 
was also known ||TTJ that departures from such an evolutionary scheme could be 



obtained with the addition of scalar fields). In fact, many of those results were quite 
interesting and contributed much to a sustained interest in the Mixmaster dynamics 
during the eighties. 

Interest in the Mixmaster universe suddenly peaked with the appearence of the 



first works on numerical experiments in 1989-90 (T^|. Those results (and others 



which followed |]T3|) were conflicting in the sense that in most cases the standard 
picture M was challenged to the effect that many workers in the field felt that a 
re-examination of the original conclusions concerning the existence of chaos in the 
Mixmaster approach to the initial singularity was necessary. At the same time other 
works appeared which either criticized [|TJ[] the use of some essentially coordinate 



dependent measures of chaoticity or indicated [15, 16] that some sort of chaotic 



behaviour should be present in the Mixmaster dynamics. 



It was felt necessary |L7j that perhaps an analytical approach, known as Painleve 
analysis, which did not share the 'defects' of numerical work could lead to more 
reliable results concerning the existence or non-existence of chaotic behaviour in the 



Mixmaster dynamics. Initial results [18] pointed to the direction of integrability 



whereas later it was realised [|19|] that the situation was more complex. It is now 
understood that this analysis cannot be used to obtain reliable results concerning 
the question of chaoticity in this model. 

Very recently Cornish and Levin |2(J] using fractal methods resolved the long- 
standing debate concerning the issue of chaoticity in the Mixmaster universe, showing 
that the system is indeed chaotic. Their analysis not only confirms the earlier ergodic 
results of Barrow || but quantifies the chaotic behaviour of the model by calculating 
a special set of numbers (topological entropy, multifractal dimension etc) relevant 
to the true dynamics. It therefore appears that all ambiquities concerning this issue 
have now disappeared. 

It is perhaps encouraging that several issues about the Mixmaster dynamics not 
connected to the question of chaos have been studied by several authors. It is indeed 



true that many problems in the dynamics of this model still remain. First of all 
it is still unknown how local the BKL analysis is pi]. Another problem is borne 



out of earlier results of Moncrief |22[ and is related to the role that the Geroch 
transformation plays for the true dymanics. Still another issue has to do with the 
description of the Mixmaster universe from the point of view of dynamical systems 
theory [p3|. We therefore believe that the rich Mixmaster behaviour will continue to 
attract the interest of cosmologists for some time to come. 

FIGURE 1. TO BE INSERTED HERE 

The purpose of this paper is to provide a rigorous analysis of the role and uses of 
adiabatic invariants in the Mixmaster problem by carefully examining the adiabatic 
invariant commonly used and introducing a new and in a sense improved invariant 
for the global Mixmaster dynamics. We then reformulate the main characteristic of 
the dynamics via a new, simpler technique along the lines of Catastrophe theory. 
We believe that further exploitation of this formulation may help to unravel certain 
global dynamical properties of the Mixmaster universe. 

The plan of this paper is as follows. In Section II we introduce our new adiabatic 
invariant for the Mixmaster system which behaves better than the standard one and 
is explicitly time-independent in the appropriate coordinates. We also perform a 
numerical simulation of the corner-run part of the evolution and we give an inter- 
pretation of the results using Misner's Hamiltonian picture. In Section III, we apply 
Catastrophe Theory as a means to gain a better understanding of the complicated 
behaviour of the model. The main result of this Section is that the passage from 
one era to the next corresponds, in the language of Catastrophe Theory, to a fold 



catastrophe which in turn may provide a potentially new way to view the global 
evolution. In the last Section, we compare our results to previous work and point 
out how generalizations to higher dimensions could be obtained. 

II Adiabatic Analysis 

Subject to a couple of overall approximations the motion of the 'universe point' 
(3 = (/3 + , (3-) is governed by the Mixmaster Lagrangian 

L full = 1A5 ((3*1 + /?f ) - 2A-k- 4n V((3), (1) 

where j3 + and /?_ are related to the shape parameters of the model, the volume 
parameter, Q plays the role of time, primes denote differentiation with respect to Q 
and A is a function of Q which evolves according to 

A' = -4e- 4a V(/?). (2) 

Here V(/3) is the standard, curvature anisotropy (Mixmaster) potential given by 

V(J3) = ^e" 4/3 + - -e"^ cosh( v / 3/5_) + ^e 2/3+ [cosh(2 v / 3/3_) - 1] + 1. (3) 

There is also the energylike equation 

4 = fi% + (3^ + 4A- 1 e~ 4f V(/3). (4) 

We begin our analysis by re-examining the role and consequences of the adiabatic 
invariant used implicitly (or explicitly) in most Mixmaster calculations. According to 
Misner the most important asymptotic form of the Mixmaster potential is when 
| /3- | is small and slowly varying and /3+ — ► +oo , that is, V(/3) ~ 1 + 4/3 2 _e 2f3+ . In 



this case Lf u u can be considered as the Lagrangian for the /?_ motion; that is, it can 
be supposed approximately that there exists a system with one degree of freedom, 
P+, depending on the slowly varying parameter /3_. This reduced system is described 
by the lagrangian, 

deduced = |A5 ^ - 8A" 2 e 2 ^" 4 "^ (5) 

in which a term 0(e _4n ) has been neglected and A and /3 + are again functions of Q 
through the solutions of @ and @ with the potential as modified. The Lagrangian 
(§) is that of a time-dependent oscillator of slowly varying frequency given by 

c^_ = 4A-!e /3+ - 2n . (6) 

Then, for the reduced system, there exists the adiabatic invariant 

E = — 

= lAe 2n -^((3'2 + u 2 _(3i), (7) 

8 

where E_ are the energy level sets and 27rE represents the area of appropriate do- 
mains bounded by curves passing through points in (the two dimensional phase) 



space (7r_,/5_). It can be shown, by adapting the methods of p4 |, that this is also 
an adiabatic invariant of the 'full' Mixmaster system (two degrees of freedom) with 
/3- slowly varying (but not necessarily small). We stress this point since it is impor- 
tant to remember that the motion of the universe point described by L reduced is only 
approximately true and the true dynamics in this case should be thought of as that 
given by Lf u u with j3- slowly varying. 

The adiabatic invariant (0) (considered originally by Misner 0) is exactly the 
one proposed by Lorentz at the first Solvay Congress in 1911 [E3, Eol. A precise 



mathematical statement concerning its range of validity was given by Littlewood 
27], in the sixties. Littlewood showed that (in our notation) 



(%)£ = c + 0(e); 

(ii) £ = c + 0(e 2 ) (£ is the average of £ over the local period 27r/u;_.) 

(Hi) there is no improvement over (i) or (ii) 

(iv) E(oo) — £(— oo) = 0(e n ) for some specific n 

provided u_ satisfied certain assumptions which give a measure, e, to the expression 
'slowly varying'. These assumptions are 

cj_ > bn, uj_ = — ■ >• as t —> ±00 (n > 1) 

dt n K ~ ' 

|c/ n) | < b n e n (n > 1), r \uo M \dt < Ci£ n_1 , (8) 

where the 6s are constants. A consequence is that, since 

r' 
\uj-{t') - w_(r)| < / |cJ_|dt -^ 

as r, t' — > —00, +00 respectively, we have 

c<j_(— 00) = c<j_(+oo). 

Strictly we need not go as far as ±00, but the in-channel time is supposed to be 
long. 

Following Arnol'd |28|] a function I(q,p;et) is an adiabatic invariant of a Hamil- 



tonian system 



dH dH , 

p=-— q=— H = H(q,p;et) (9) 



if Vk > 3 e such that if e < e and < t < l/e , 

|/(g(t),p(t);et)-/(g(0),p(0);0)|<«. (10) 

The action variable of the corresponding autonomous problem is always an adiabatic 
invariant and it is for this reason that (0) can be selected as the adiabatic invariant. 



However, the perpetuality (cf. f24fl ) must be calculated and verified. Previous 
works provide no information about the behaviour of uj_ and certainly no measure 
of e. Fortunately it is a fairly straightforward matter to test the validity of the 
use of the adiabatic invariant (0). We can simply integrate the equations of motion 
numerically and substitute the numbers into (0) to observe the variation of S with 

n. 

In the channel regime the equations governing the motion of the system point are 

A' = _i6e^+- 4C /?! (11) 

(3' + = ±(4-(^) 2 -^- 4 ^ 2 ) 1/2 (12) 

/T - |e^+- 4 W + I^+- 4f V?_ = 0. (13) 

A A 

The positive sign in ( fl2j) applies until the expression under the square root sign be- 
comes zero. Equations (|TTD-(|T3|) are those of Misner except that (3q has not been 
introduced. We take initial conditions consistent with the assumptions governing 
motion in the channel. We integrate the system (|TT|-|i"3"D numerically using an imple- 



mentation of the Runge-Kutta scheme p9|. The results are illustrated in Figure 2 



It is quite evident that these computations do not support the use of the adiabatic 
invariant over the whole time spent inside the channel. 
FIGURE 2. TO BE INSERTED HERE. 



Let us now consider some consequences of the above analysis. Upon the introduc- 
tion of a relative coordinate 

A, = /?+ - 2Q (14) 

and, with the help of eqs © and (0), eq (§) becomes 

= /% + 4& + ?ffe-*>. (15) 

On the assumption that (3' is small, (3® is neglected and (p~5|) gives immediately 

fa = log(fi - 0) + const (16) 

which leads to /3o — > — oo as fi approaches the critical value Q,q. This means 
that the particle leaves the channel and returns bouncing in the triangular region. 

Consider now the assumption that is usually made namely, that (3'1 may be ne- 
glected. We claim that this assumption leads to a solution which is only asymptot- 
ically correct. To see this, notice that it follows easily that the solution of ( |I5|) is 
given implicitly by 

n - n = - la~ l e\ Po + y/a- 2 eP° - l| - \ log {a~ x A^ + y/a~ 2 e^ - l| , (17) 



where 



.2 



(18) 



is assumed constant. Misner's original solution (|16"D is, in a sense, asymptotic to (|T7D 
and our assertion follows. 

Of more interest, however, is the actual equation for /3' Q , which follows immediately 
from (|T5D, viz. 



^ = -2±2 v / l-a 2 e-/ 3 o. (19) 



8 



(The upper sign applies for the initial motion in the channel since \(3' + \ < 2.) Eq. 
(O) is valid only if 

e?° > a 2 . (20) 

Inserting this into eq ([[7]) we find that 

Q<Q -- or (3 > -0.60 + const. (21) 

In other words Q never reaches the critical value Qq and /5o cannot go to — oo which 
is what is necessary for the particle to return to the triangle from the channel. (Of 
course, in practice, it is only sufficient to have /?o large and negative for the particle to 
resume bouncing with the wall in the triangular box as soon as it leaves the channel, 



but even for this to happen the bound in Eq. ([ZT]) seems too stringent.) 

We are now ready to introduce a new invariant which in many ways behaves better 
than S. It is a well-known fact that the time-dependent oscillator described by the 
Hamiltonian 



H = W + u\t)q 2 ) (22) 



possesses the first integral 



I=\\(PP-PQ) 2 +(j)) ), (23) 

which is known as the Ermakov-Lewis invariant |30]| , provided that the auxiliary 

variable is a solution of 

p + u 2 (t)p=^. (24) 

P 3 



Furthermore the solution of (0) has been given by Pinney |3T| in terms of the linearly 
independent solutions of 

v + LU 2 (t)v = 0. (25) 



9 



The Lagrangian (|5]) gives directly the Hamiltonian 

deduced = |A^n 2 _ + 8A^e 2/3 +- 4C /3!, (26) 

where 

n = A^P'_. (27) 



The reduced Hamiltonian fl26|) is not precisely of the form of (p2|) since the coefficient 
of Tl 2 _ is not constant and we have the equivalent of an harmonic oscillator of variable 
mass. However, following the procedure detailed by Leach [02], for the treatment of 



a Hamiltonian of the form of (|26"D we introduce a change of time scale 

T= f A~2(fT)dfi (2* 

so that the Mixmaster system in this regime is described by 



H reduced = |n 2 _ + \u 2 {T)(3 2 _ (29) 



in which T is now the independent variable and 

cj 2 (T(fi)) = i6e 2/3 +( Q )- 4Q . (30) 

Under the generalised canonical transformation [33 

Q=^— P = pp-pq T=[p- 2 (T)dT (31) 

P J 

(^9|) is transformed to 

H=\{P 2 + Q 2 ) (32) 

provided that p{T) is a solution of (|2~4D (with cu 2 (T) from (|30D instead of the ui 2 {t)). 



H is the Ermakov-Lewis invariant in the canonical variables in which it becomes free 
of explicit dependence on time. 

10 



It is H which should be used instead of the adiabatic invariant used by Misner. 
However, there is a problem. In the new variables the evolution of the oscillator 
is easily described, but that of j3 + becomes unmaneageable as the Ermakov-Lewis 
invariant does not lead to a significant simplification of ([12]). We must resort to 
numerical computation and this may as well be performed in the original variables. 

FIGURES 3,4,5. TO BE INSERTED HERE 

Over the interval < uj < 30000, (3 + increases essentially linearly with u. There is 
no indication that it approaches — oo and as the deviation from strict linearity is so 
small, it can only be expected to take a long time to approach zero. In the meantime 
the potential well proceeds outwards as is evident from Figure 4 at which the contour 
is that corresponding to the energy of the motion. We see that the wall outmarches 
the (3+ value of the universe point. Over this period, as depicted in Figure 5, the 
amplitude of the (3- motion increases in a strictly monotonic fashion with a rate of 
increase increasing with time. 

It is clear from these results that, as the particle moves along the channel, the 
potential walls are receding. Initially this does not present a problem as the /?+ 
velocity of the particle is sufficiently large. However, there is a critical time f2 when 
the walls 'leave the particle behind'. Hence the particle finds itself in the triangular 
region again and no longer in the channel. Note that, as f3 + is never less than zero, 
the particle does not 'turn around' in the sense that is sometimes described. 



11 



Ill A Catastrophe Description 

In what follows we present a novel way to describe the qualitative differences of 
the main stages in the evolution of the Mixmaster universe. As we shall see, the 
standard interpretation can be reached quite naturally and independently by this line 
of thought. This approach is established by the application of Singularity Theory 



fl34jl and in particular by that branch of singularity theory known as Catastrophe 
Theory |J3|. Catastrophe theory studies changes in the equilibria of potentials as 
the control parameters of the system change. The local properties of the potential 
in a gradient or a dynamical system are determined by a sequence of theorems such 
as the Implicit Function Theorem of advanced calculus, the Morse Lemma |5^] and 



the Thorn Theorem [37] 



We consider the Mixmaster universe as a gradient system described by the poten- 
tial (^). We choose as a control parameter the volume (time) parameter Q. We set 
V = (d/dfl + , d/df3_) and denote the Hessian matrix by 



Va 



( d 2 V d 2 V \ 

d 2 V d 2 V 



(33) 



Then, for the Kasner-to-Kasner evolution described by the potential 

P+ — -oo, V(p) ~ \e~^\ (34) 

we find 

VV=(--e- 8l3+ ,0) ^0. (35) 

This means that during an era the Implicit Function Theorem applies and there are 
no critical points. (This also implies that there is a smooth change of coordinates 



12 



which makes the potential (^|) depend on only one of the variables, say j3 + . We see 
that the form of the potential (B3) can be deduced from the general form (^) by 
using this argument without resorting to any sort of approximations.) 

Secondly we examine the structure of the potential in the neighborhood of the 
isotropy point (0,0) given by 

(/?+,/?_) ~ (0,0), V{(3) ~ 16(# + /£), (36) 

After some straightforward manipulations we find that 

W = 0. (37) 

We see that the conditions for the validity of the Implicit Function Theorem are 
no longer satisfied. Equation ( |3"TD implies that near the isotropic point (0, 0) the 
Mixmaster universe is in a stable equilibrium state. To see this we determine the 
stability properties of this state by finding the eigenvalues of the Hessian matrix, Vij. 
Firstly after a tedious calculation we find 

256 
det m Vij = — + 0. (38) 

Also the eigenvalues of V^ are Ai = A2 = 16. This means that due to the Morse 



theorem [|36] there is a smooth change of variables so that the potential in this case 
takes the form 

V = Ml = \\ + \\ = 16{f3 2 + + (3 2 _). (39) 

Not surprisingly this is exactly the form of the potential that Misner found in this 
case. Mq stands for the Morse 0-saddle which is the only i-saddle that is stable 
for two-dimensional gradient systems (c/ [0). Thus the point (0,0) is a Morse 
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critical point (isolated, nondegenerate) . In particular, the potential in this case is 
structurally stable. 

Lastly we examine the corner-run evolution which turns out to be the most in- 
teresting from the point of view of Catastrophe Theory. In this case we find 

W = (40) 

and 

detVij = -8172/?!e 8/3+ . (41) 

It is clear that for all points on the /3 + -axis (/3_ = 0) we have 

det^^Vij = 0. (42) 

This implies that in the channel region all points which lie on the /3 + -axis are non- 
Morse critical points. In this case we can cast the potential in a canonical form 



by adopting a procedure known as the Thorn Splitting Lemma [p8| . We split the 
potential into a Morse part and a non-Morse part according to the number of the 
vanishing eigenvalues of the Hessian matrix for this case. These are found to be 

Ai = 0, A 2 = 32e 4/3+ , (43) 



which due to a theorem of Thorn [37j guarantees that there is a smooth change of 



variables that puts the potential (in the channel) in the decomposed form 

V((3)=(3 3 + + ap + + 32(3 2 _, (44) 

with a^O. The first two terms in this potential (the non-Morse part) form what 
is known as the fold Catastrophe (A 2 ) and it is the simplest of the seven elemen- 
tary catastrophes first discussed by Thorn in p7| . The Morse part of the above 

14 



decomposition is unaffected by perturbations so it is only necessary to study how 
the qualitative properties of the catastrophe function A 2 = f3\ + a[3 + are changed 
as the control parameter changes. When a > there are no critical points whereas 
a < gives two critical points namely, (3± = ±y/—a. The case a = is the separa- 
trix in the control parameter space between functions of two qualitatively different 
types (no critical points and two critical points). 

Our interpretation of the above results uses the delay convention of Catastrophe 



Theory (see for instance [35]). Since, as we have shown, during the bounces of the 
point with the walls in the triangular box there are no critical points, we imagine that 
when the point enters the channel has a < 0. Then as it moves inside the channel 
the degenerate point a = is reached (the stable minimum disappears into the 
degenerate critical point). At this instant /?_ is no longer small, there are no critical 
points and the system jumps to the lowest of the two minima (the stable attractor) 
of a > 0. This produces a (point) shock wave which is the simplest elementary 
catastrophe (fold). This, in turn, means that the system (point) has found itself 
bouncing again inside the (now larger) triangular box. 

IV Conclusions 

Our adiabatic analysis relates to the well-known issue of the so-called 'anomalous' 
behaviour discussed previously analytically by Berger in Ref. [13|] and numerically 



in ref. |39|] , [§tj] . Physically, this in-channel behaviour appears only when the initial 
value of the so-called BKL parameter u is sufficiently large (a BKL 'long era'). 
However, the required value of u becomes larger (and therefore less probable) as 
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the singularity is approached. This in turn means that the in-channel behaviour 
becomes less probable as the Mixmaster singularity is approached. 



The recent demonstration of chaoticity by Cornish and Levin |20| via the existence 
of a Mixmaster fractal strange repellor may be seen in the light of the non-adiabatic 
Mixmaster evolution discussed here. It is interesting to point out that the problem 
of the existence of an adiabatic invariant for higher dimensional generalizations of 
the Mixmaster universe II or in higher derivate extensions JT0J (wherein chaotic 



behaviour may be absent) is a nontrivial one and one expect that the usual difficulties 



24] present in dynamical systems with more than two degrees of freedom exist in 



this problem too. 



Our numerical results parallel those given in [39, [5]] in the following respects: In 
those references, figures equivalent to our Fig. 5 are given but the variables (3±/fl 
are plotted there (mixing bounces) rather than our variables j3±. We stress that 
no confusion must arise in this respect since, in the former variables the trajectory 
associated with a single era appears to move outward along a corner and then inward 
again while, here the motion is strictly outward. Further, as is clearly emphasized 



by Berger in pO ], the angle of the minisuperspace trajectory becomes ever closer 
to the perpendicular to the ray down the corner as the era progresses towards the 
singularity. A change of era occurs when the trajectory points inwards rather than 
outwards with respect to this perpenticular direction. 

We hope that our reformulation of the problem in terms of Catastrophe theory in 
Section III may be further used to examine questions of current interest such as, for 
instance, issues connected with the occurence of chaotic behaviour. In some sense, 
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our Catastrophe results correspond to just Mixmaster statics. Further dynamical 
issues could be addressed if one considers the Mixmaster system as a gradient system 
as is usual in Catastrophe discussions of dynamical systems. 

Another issue that is raised by our formulation is the influence of changing the 
time parameter on the character of the Catastrophe. It is well-known that there exist 
different time parametrizations for the description of the Mixmaster dynamics (see, 



for instance, |4T)|). It is therefore appropriate to ask how the Catastrophe profile 
of the Mixmaster dynamics is affected by different choices of time. Although the 
answer to this question is uncertain at present, we believe that a physically relevant 
formulation should be unaffected by different time choices. 
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